using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.QuadPack
{
    /// <summary>
    /// DQAGE - Approximation to definite integral. Allows user's choice of Gauss-Kronrod integration rule.
    /// </summary>
    [Serializable]
    public static class DqageClass
    {
        /// <summary>
        /// DQAGE - Approximation to definite integral. Allows user's choice of Gauss-Kronrod integration rule.
        /// </summary>
        /// <param name="f">Double precision function to be integrated.</param>
        /// <param name="a">Lower limit of integration.</param>
        /// <param name="b">Upper limit of integration.</param>
        /// <param name="epsabs">The absolute accuracy requested.</param>
        /// <param name="epsrel">The relative accuracy requested.</param>
        /// <param name="irule">The integration rule to be used as follows: irule = 1 -- G_K 7-15; irule = 2 -- G_K 10-21; irule = 3 -- G_K 15-31; irule = 4 -- G_K 20-41; irule = 5 -- G_K 25-51; irule = 6 -- G_K 30-61.</param>
        /// <param name="abserr">The requested absolute error of the result.</param>
        /// <param name="neval">The needed function evalations.</param>
        /// <param name="ier">The requested errorcode.</param>
        /// <param name="last">The maximum number of subintervals.</param>
        /// <returns>The result of the numerical integration.</returns>
        public static double Dqage(HardRealFunction f, double a, double b, double epsabs, double epsrel, int irule,
                                   ref double abserr, ref int neval, ref int ier, ref int last)
        {
            double area;
            double area1 = 0;
            double area2 = 0;
            double area12;
            double a1;
            double a2;
            double b1;
            double b2;
            double c;
            double defabs = 0;
            double defab1 = 0;
            double defab2 = 0;
            double errbnd = 0;
            double errmax = 0;
            double error1 = 0;
            double error2 = 0;
            double erro12 = 0;
            double errsum = 0;
            double resabs = 0;
            double result = 0;
            double[] alist = new double[QuadPackConstants.LIMIT];
            double[] blist = new double[QuadPackConstants.LIMIT];
            double[] rlist = new double[QuadPackConstants.LIMIT];
            double[] elist = new double[QuadPackConstants.LIMIT];
            int iroff1;
            int iroff2;
            int k;
            int keyf;
            int maxerr;
            int nrmax;
            int[] iord = new int[QuadPackConstants.LIMIT];
            int limit;

            limit = QuadPackConstants.LIMIT - 1;
            ier = 0;
            neval = 0;
            last = 0;
            result = 0.0;
            abserr = 0.0;
            alist[0] = a;
            blist[0] = b;
            rlist[0] = 0.0;
            elist[0] = 0.0;
            iord[0] = 0;
            defabs = 0.0;
            resabs = 0.0;
            if ((epsabs < 0.0) && (epsrel < 0.0))
            {
                ier = 6;
            }
            if (ier == 6)
            {
                return result;
            }

            // First approximation to the integral. 
            keyf = irule;
            if (irule <= 0)
            {
                keyf = 1;
            }
            if (irule >= 7)
            {
                keyf = 6;
            }
            c = keyf;
            neval = 0;
            switch (keyf)
            {
                case 1:
                    result = Dqk15Class.G_K15(f, a, b, ref abserr, ref defabs, ref resabs);
                    break;
                case 2:
                    result = Dqk21Class.G_K21(f, a, b, ref abserr, ref defabs, ref resabs);
                    break;
                case 3:
                    result = Dqk31Class.G_K31(f, a, b, ref abserr, ref defabs, ref resabs);
                    break;
                case 4:
                    result = Dqk41Class.G_K41(f, a, b, ref abserr, ref defabs, ref resabs);
                    break;
                case 5:
                    result = Dqk51Class.G_K51(f, a, b, ref abserr, ref defabs, ref resabs);
                    break;
                case 6:
                    result = Dqk61Class.G_K61(f, a, b, ref abserr, ref defabs, ref resabs);
                    break;
            }
            last = 0;
            rlist[0] = result;
            elist[0] = abserr;
            iord[0] = 0;

            // Test on accuracy. 
            errbnd = (((epsabs) > (epsrel * Math.Abs(result))) ? (epsabs) : (epsrel * Math.Abs(result)));
            if ((abserr <= 50.0 * Double.Epsilon * defabs) && (abserr > errbnd))
            {
                ier = 2;
            }
            if (limit == 0)
            {
                ier = 1;
            }
            if ((ier != 0) || ((abserr <= errbnd) && (abserr != resabs)) || (abserr == 0.0))
            {
                goto _60;
            }

            // Initialization. 
            errmax = abserr;
            maxerr = 0;
            area = result;
            errsum = abserr;
            nrmax = 0;
            iroff1 = 0;
            iroff2 = 0;

            // Main Loop. 
            for (last = 1; last <= limit; (last)++)
            {
                // Bisect the subinterval with the largest error estimate. 
                a1 = alist[maxerr];
                b1 = 0.5 * (alist[maxerr] + blist[maxerr]);
                a2 = b1;
                b2 = blist[maxerr];
                switch (keyf)
                {
                    case 1:
                        area1 = Dqk15Class.G_K15(f, a1, b1, ref error1, ref resabs, ref defab1);
                        area2 = Dqk15Class.G_K15(f, a2, b2, ref error2, ref resabs, ref defab2);
                        break;
                    case 2:
                        area1 = Dqk21Class.G_K21(f, a1, b1, ref error1, ref resabs, ref defab1);
                        area2 = Dqk21Class.G_K21(f, a2, b2, ref error2, ref resabs, ref defab2);
                        break;
                    case 3:
                        area1 = Dqk31Class.G_K31(f, a1, b1, ref error1, ref resabs, ref defab1);
                        area2 = Dqk31Class.G_K31(f, a2, b2, ref error2, ref resabs, ref defab2);
                        break;
                    case 4:
                        area1 = Dqk41Class.G_K41(f, a1, b1, ref error1, ref resabs, ref defab1);
                        area2 = Dqk41Class.G_K41(f, a2, b2, ref error2, ref resabs, ref defab2);
                        break;
                    case 5:
                        area1 = Dqk51Class.G_K51(f, a1, b1, ref error1, ref resabs, ref defab1);
                        area2 = Dqk51Class.G_K51(f, a2, b2, ref error2, ref resabs, ref defab2);
                        break;
                    case 6:
                        area1 = Dqk61Class.G_K61(f, a1, b1, ref error1, ref resabs, ref defab1);
                        area2 = Dqk61Class.G_K61(f, a2, b2, ref error2, ref resabs, ref defab2);
                        break;
                }

                // Improve previous approximations to integral and error,
                //		and test for accuracy. 
                (neval) += 1;
                area12 = area1 + area2;
                erro12 = error1 + error2;
                errsum = errsum + erro12 - errmax;
                area = area + area12 - rlist[maxerr];
                if ((defab1 != error1) && (defab2 != error2))
                {
                    if ((Math.Abs(rlist[maxerr] - area12) <= 1.0e-5 * Math.Abs(area12)) && (erro12 >= .99 * errmax))
                    {
                        iroff1++;
                    }
                    if ((last > 9) && (erro12 > errmax))
                    {
                        iroff2++;
                    }
                }
                rlist[maxerr] = area1;
                rlist[last] = area2;
                errbnd = (((epsabs) > (epsrel * Math.Abs(area))) ? (epsabs) : (epsrel * Math.Abs(area)));
                if (errsum > errbnd)
                {
                    // Test for roundoff error and eventually set error flag. 
                    if ((iroff1 > 6) || (iroff2 > 20))
                    {
                        ier = 2;
                    }

                    // Set error flag in the case that the number of subintervals
                    //	equals the limit. 
                    if (last == limit)
                    {
                        ier = 1;
                    }

                    // Set error flag in the case of bad integrand behavior at a
                    //	point of the integration range. 
                    if ((((Math.Abs(a1)) > (Math.Abs(b2))) ? (Math.Abs(a1)) : (Math.Abs(b2))) <=
                        (1.0 + c * 1000.0 * Double.Epsilon) * (Math.Abs(a2) + 1.0e4 * Double.MinValue))
                    {
                        ier = 3;
                    }
                }
                // Append the newly-created intervals to the list. 

                if (error2 <= error1)
                {
                    alist[last] = a2;
                    blist[maxerr] = b1;
                    blist[last] = b2;
                    elist[maxerr] = error1;
                    elist[last] = error2;
                }
                else
                {
                    alist[maxerr] = a2;
                    alist[last] = a1;
                    blist[last] = b1;
                    rlist[maxerr] = area2;
                    rlist[last] = area1;
                    elist[maxerr] = error2;
                    elist[last] = error1;
                }

                // Call DQSORT to maintain the descending ordering in the list of
                //	error estimates and select the subinterval with the
                //	largest error estimate (to be bisected next). 

                DqsortClass.Dqsort(limit, last, ref maxerr, ref errmax, elist, iord, ref nrmax);
                if ((ier != 0) || (errsum <= errbnd))
                {
                    break;
                }
            }

            // Compute final result. 

            result = 0.0;
            for (k = 0; k <= last; k++)
            {
                result += rlist[k];
            }
            abserr = errsum;
            _60:
            if (keyf != 1)
            {
                neval = (10 * keyf + 1) * (2 * (neval) + 1);
            }
            else
            {
                neval = 30 * (neval) + 15;
            }

            return result;
        }
    }
}